function r = sh2(N)

close all
  

  
  %time_steps = [0:N-1]./(N-1)*1;
  time_steps = chebpoints(N, 0, 1);

  function resB = B(z)
    if z == 0
      resB = 6;
    else
      resB = (1/3).*z.^(-1).*(1+z+z.^2).^(-2).*(2.*3.^(1/2).*pi+6.*z+3.^(1/2).*pi.*z+6.*z.^2+3.*z.^3+(-2).*3.^(1/2).*(2+z).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*(2+z).*log(1+z+z.^2));
    end
  end
  function resSk = Sk(z)
    if (z < 0.001)
      resSk = (14/3)+(25/3).*z+(-22/3).*z.^2+(-7/3).*z.^3+(221/21).*z.^4+(-325/42).*z.^5+(-244/63).*z.^6+(7741/630).*z.^7+(-5167/630).*z.^8+(-217/45).*z.^9+(12323/910).*z.^10;
    else
      %resSk = z.^3.*(1+z+z.^2).^(-1).*((-1)+4.*z.^(-3)+4.*z.^(-2)+2.*z.^(-1)+(1/3).*z.^(-4).*((-4)+z.*(1+z).*((-6)+z.*(2+z.*(2+z)))).*(3.^(1/2).*pi+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2)));
      %resSk = z.^(-1).*(4+z.^3.*(1+z+z.^2).^(-1).*((-1)+4.*z.^(-3)+4.*z.^(-2)+2.*z.^(-1)+(1/3).*z.^(-4).*((-4)+z.*(1+z).*((-6)+z.*(2+z.*(2+z)))).*(3.^(1/2).*pi+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2)))) + 4;
      resSk = z.^(-1).*(4+z.^(-1).*(4+z.^3.*(1+z+z.^2).^(-1).*((-1)+4.*z.^(-3)+4.*z.^(-2)+2.*z.^(-1)+(1/3).*z.^(-4).*((-4)+z.*(1+z).*((-6)+z.*(2+z.*(2+z)))).*(3.^(1/2).*pi+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2)))));

    end
  end
  function resM = Mhk(z)
    resM = zeros(3,3);
    resM(1,1) = 1;
    resM(2,1) = 4;
    resM(2,2) = z;
    %resM(3,1) = 4-z.^3; %uncomment if solving mixed pde
    resM(3,3) = 2;
  end
  function res = deqhk(z, u)
    h   = u(1);
    hz  = u(2);
    k   = u(3);

    res = zeros(2,1);
    %kz = Sk(z) + z + 3.*z.^2.*h; % uncomment if solving mixed pde
    kz = Sk(z);% comment out if solving mixed pde

    
    res(1) = hz;
    res(2) = B(z);
    res(3) = kz;
  end
  uin = [0; 6./4; 0];
  opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10, 'MStateDependence', 'none', 'Mass', @Mhk);
  [tout, uout] = ode15s(@deqhk, time_steps, uin, opts);
  z = time_steps';
  nh2b = uout(:, 1);
  nh2b = 1/2 + z.*nh2b; % with z^2 factored out
  nk2d = uout(:, 3);
  %nk2d = 1 + z.*nk2d; % uncomment if solving mixed pde
  nk2d = 2 + z.*nk2d;  % comment out if solving mixed pde
  nk2d = nk2d - z.*log(z)/2; % with 1/z^2 factored out

  save('nf.m', 'nh2b', 'nk2d');

  plot(z, nh2b, z, nk2d);

  % tests 


end
